from opnexrad import opnexrad
from opaod import *
from opaodbase import opaod_base
from opecmwf import *
from oppopdens import *
from oplandcover import *
from opsoil import *
from opglim import *
from oplith4 import *
from opgebco import *
from opsa import *
from opmodel import *
from mpl_toolkits.axes_grid1 import make_axes_locatable
import xarray as xr
import os
from cartopy.feature import ShapelyFeature
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import pyart
import warnings
warnings.simplefilter('ignore')
import time
import importlib
import joblib
%load_ext autoreload
%autoreload 2
The autoreload extension is already loaded. To reload it, use: %reload_ext autoreload
base = "/home/xiaohe/Documents"
fname = os.path.join(base,'NEXRAD/AirNow/US_Boundary/USBoundary.shp')
shape_feature = ShapelyFeature(Reader(fname).geometries(),
ccrs.PlateCarree(), facecolor='none')
AOD_UTC = "2020-01-06T18:30"
mdl_type ="nonex" #"nex" , "nonex" both nexrad and ecmwf,
if mdl_type == "nonex":
print("Load Mdl")
Mdl = joblib.load(os.path.join(base,"NEXRAD/Model/ResidualAnalysis/US_AODDQF_SA_Z01_ET_Bayes_testMdl/ETreg_ECMWF_random.joblib"), mmap_mode=None)
predictors = [
'uv10', 'd2m', 't2m', 'lai_hv', 'lai_lv', 'sp', 'sro', 'tp',
'popden', 'landcover', 'soil', 'glim', 'lith4', 'gebco','AOD','DQF','LZA','SZA','SAA']
### AOD base ################################################################################
ds_GOES = xr.open_dataset("noaa-goes16/ABI-L2-AODC/2020/006/18/OR_ABI-L2-AODC-M6_G16_s20200061831177_e20200061833550_c20200061840213.nc",decode_times=False)
print("trim and grid AOD, save to df")
df_NG = opaod_base(ds_GOES,plotvar="yes",shape_feature=shape_feature)
full_extent = df_NG.iloc[:,0:2] # used as base for predicted model merge
filter_NG = df_NG.dropna()
elif mdl_type == "nex":
print("Load Mdl")
Mdl = joblib.load(os.path.join(base,"NEXRAD/Model/ResidualAnalysis/US_AODDQF_SA_Z01_ET_Bayes_testMdl/ETreg_random.joblib"), mmap_mode=None)
predictors = [
'z01reflectivity', 'z01velocity',
'z01spectrum_width', 'z01differential_phase',
'z01differential_reflectivity', 'z01cross_correlation_ratio',
'uv10', 'd2m', 't2m', 'lai_hv', 'lai_lv', 'sp', 'sro', 'tp',
'popden', 'landcover', 'soil', 'glim', 'lith4', 'gebco','AOD','DQF','LZA','SZA','SAA']
## Read NEXRAD file and retrived iregular 2d coordinates
grid = pyart.io.read_grid(os.path.join(base,"NEXRAD/Variables/TestingNEXRAD/2020-01-06_183000.nc"))
NEXZ01, full_extent = opnexrad(grid,shape_feature,plotvar="reflectivity")
## AOD match ################################################################################
stamp = "2020-01-06T18:00"
filter_N = NEXZ01.dropna() # if nexrad is included
start_time = time.time()
nc_path = match_goes(filter_N, stamp)
filter_NG = filter_N.dropna()
print("-------%s seconds -------" %(time.time() - start_time))
print("potting AOD")
aod_plot(nc_path,shape_feature,"AOD")
Load Mdl trim and grid AOD, save to df -------3.1011202335357666 seconds -------
### ECMWF #######################################################################################
ds = xr.load_dataset(os.path.join(base,'NEXRAD/Model/Operation/era5_US_Land_2020_1_6_1800.grib'), engine='cfgrib')
match_ecmwf(filter_NG,ECMWF_xr=ds)
filter_NEG = filter_NG.dropna()
# plot emcwf
plot_2decmwf(ds,shape_feature)
retrieving u10 values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving v10 values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving d2m values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving t2m values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving lai_hv values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving lai_lv values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving sp values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving sro values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF retrieving tp values from ECMWF grid file ...... There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
### popdens ####################################################################################
start_time = time.time()
usda = combinetiff("PopDens/us_pop15.tif",
"PopDens/us_pop20.tif")
# only support pop from 15 to 23
match_popdens(filter_NEG,POPDEN_xr=usda, year="2020")
print("-------%s seconds -------" %(time.time() - start_time))
# ## plot pop dens tif
da20 = xr.open_rasterio("PopDens/us_pop20.tif")
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da20.variable.data[0],vmin=0,vmax=100)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("Popden2020.png")
filter_NEG_P = filter_NEG.dropna()
retrieving population density values in 2020 from CEDAC geotiff file ...... There are 20849 sites macthed with the population density raster, 0 sites are out of coverage from Radar -------36.47853183746338 seconds -------
#landcover##########################################################################################
start_time = time.time()
da = xr.open_rasterio("LandCover/Reclass_NLCD_Re1kmMajorit.tif")
match_landcover(filter_NEG_P,da)
print("-------%s seconds -------" %(time.time() - start_time))
## plot
lc = xr.open_rasterio("LandCover/Reclass_NLCD_Re1kmMajorit.tif")
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(lc.variable.data[0],vmin=1,vmax=18)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("LandCover.png")
filter_NEG_PL = filter_NEG_P.dropna()
### soil ##########################################################################################
start_time = time.time()
da = xr.open_rasterio("Soil/so2015v2.tif")
print("Triming to US extent")
da = da.where(da.y<50,drop=True).where(da.y>25,drop=True).where(da.x>-125,drop=True).where(da.x<-65,drop=True)
match_soil(filter_NEG_PL,da)
print("-------%s seconds -------" %(time.time() - start_time))
filter_NEG_PLS = filter_NEG_PL.dropna()
## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=1,vmax=100)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("soil.png")
### glim ##########################################################################################
start_time = time.time()
da = xr.open_rasterio("GLiM/glim_wgs84_0point5deg.tif")
print("Triming to US extent")
da = da.where(da.y<50,drop=True).where(da.y>25,drop=True).where(da.x>-125,drop=True).where(da.x<-65,drop=True)
match_glim(filter_NEG_PLS,da)
print("-------%s seconds -------" %(time.time() - start_time))
filter_NEG_PLSG = filter_NEG_PLS.dropna()
## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=1,vmax=13)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("glim.png")
### lith4 ##########################################################################################
star_time = time.time()
da = xr.open_rasterio("Lithology4Class/fullareagmnarocktypedecde.tif")
match_lith4(filter_NEG_PLSG, da)
print("-------%s seconds -------" %(time.time() - start_time))
filter_NEG_PLSGL = filter_NEG_PLSG.dropna()
## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=1,vmax=4)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("lith4.png")
### gebco ##########################################################################################
start_time = time.time()
da = xr.open_rasterio("GEBCO_2020_tif/gebco_2020_n50.0_s25.0_w-125.0_e-65.0.tif")
match_gebco(filter_NEG_PLSGL, da)
print("-------%s seconds -------" %(time.time() - start_time))
filter_NEG_PLSGLG = filter_NEG_PLSGL.dropna()
## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=-100,vmax=3300)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("gebco.png")
retrieving landcover values from geotiff file ...... There are 20849 sites macthed with the landcover raster, 0 sites are out of coverage from Radar -------48.6298770904541 seconds ------- Triming to US extent retrieving soild values from geotiff file ...... There are 20849 sites macthed with the landcover raster, 0 sites are out of coverage from Radar -------40.850067138671875 seconds ------- Triming to US extent retrieving soild values from geotiff file ...... There are 20849 sites macthed with the landcover raster, 0 sites are out of coverage from Radar -------40.30678725242615 seconds ------- retrieving values from geotiff file ...... There are 20849 sites matched with the raster, 0 sites are out of coverage -------90.83371686935425 seconds ------- retrieving values from geotiff file ...... There are 20849 sites matched with the raster, 0 sites are out of coverage -------54.32755208015442 seconds -------
### SA ##########################################################################################
match_sa(filter_NEG_PLSGLG,UTC = AOD_UTC)
filter_NEG_PLSGLGS = filter_NEG_PLSGLG.dropna()
# plot_sa(filter_NEG_PLSGLGS,shape_feature,plotvar="SZA")
### Model ##########################################################################################
modelPredict(filter_NEG_PLSGLGS, full_extent, predictors, Mdl,shape_feature)
Calculating LZA Calculating SA Progress: [----------------------------->] 100 %
| Latitude | Longitude | predictions | |
|---|---|---|---|
| 0 | 25.0 | -125.000000 | NaN |
| 1 | 25.0 | -124.874739 | NaN |
| 2 | 25.0 | -124.749478 | NaN |
| 3 | 25.0 | -124.624217 | NaN |
| 4 | 25.0 | -124.498956 | NaN |
| ... | ... | ... | ... |
| 158395 | 50.0 | -65.501044 | NaN |
| 158396 | 50.0 | -65.375783 | NaN |
| 158397 | 50.0 | -65.250522 | NaN |
| 158398 | 50.0 | -65.125261 | NaN |
| 158399 | 50.0 | -65.000000 | NaN |
158400 rows × 3 columns
# # use 2km resoluation as base grid, need to regrid popdens landcover lith4 gebco to 2km grid
# da20 = xr.open_rasterio("PopDens/us_pop20.tif")
# ds = da20
# xx, yy = np.meshgrid(ds.x.data, ds.y.data)
# x_flat = xx.flatten()
# y_flat = yy.flatten()
# x_flat[x_flat<-125] = None
# x_flat[x_flat>-65] = None
# y_flat[y_flat>50] = None
# y_flat[y_flat<25] = None
# arr = ds.data
# arr_flat = arr.flatten()
# trimmed = pd.DataFrame()
# trimmed["Latitude"] = y_flat
# trimmed["Longitude"] = x_flat
# trimmed["popdens"] = arr_flat
# ## filter out None values which are outside US
# trimmed_drop = trimmed.dropna(subset=["Latitude"])
# trimmed_drop = trimmed_drop.dropna(subset=["Longitude"])
# ## 10km grid which is the same as nexrad extent
# xi = np.linspace(-125,-65,2400)
# yi = np.linspace(25,50,1650)
# xi,yi = np.meshgrid(xi,yi)
# zi = griddata((trimmed_drop["Longitude"], trimmed_drop["Latitude"]),trimmed_drop["popdens"], (xi,yi), method='linear')
# ds_NEXREG = xr.DataArray(zi, coords=[yi[:,1], xi[0]], dims=["lat", "lon"])
# fig = plt.figure(figsize=[20,13])
# ax = plt.subplot(projection=ccrs.PlateCarree())
# ax.add_feature(shape_feature,edgecolor='blue')
# ax.gridlines(draw_labels=True)
# ax.coastlines()
# ds_NEXREG.plot(cmap=plt.cm.Reds, vmin=0,vmax=100, x='lon', y='lat')
# plt.savefig("poptestxrlinear.png")